%% segement all the files
FileList = {
   
  %  'mohEEG111.mat'
    'Marion12_2504.mat'
  %  'Tim12_2804.mat'
   %   'Romain12_2904.mat'
    };

%FileList = { 'session2.acq.mat' };
%Fs = 2000; default
Fs = 2000;

window = 3000;
X = []; % contains the trials from all classes
Y = []; % contains the classification from all classes

is_pulsar_single_trigger_channel_from_eeglab = false;

for curFileIx = 1:length(FileList)
    
        %Instead of file the data is obtained from EEGLAB 
        if (is_pulsar_single_trigger_channel_from_eeglab)
            
            FileList = {'EEGLAB data'};
            data = EEG.data'; 
            Fs = EEG.srate;
            
            datar = data(:,1:size(data,2)-1); % remove 1 trigger channel at the end

            trigger_channel = data(:, size(data,2));

            classN = length(count_sequence(trigger_channel,2)); % "2" tags each class

            classRest = classN; % the rest class is always the last

            trialsN = length(count_sequence(trigger_channel,1));  % "1" tags each trial
            
        else   
            
            load(['D:\temp\subjects\' FileList{curFileIx}]);
            disp(['Subject file:' FileList{curFileIx}]);
            
            datar = data(:,1:size(data,2)-2); % remove 2 trigger channel at the end

            %trial_channel = data(:,8);
            trial_channel = data(:, size(data,2)-1);

            %class_separation_channel = data(:,9);
            class_separation_channel = data(:, size(data,2));

            classN = length(count_sequence(class_separation_channel,5));  % "5" tags each class
            

            classRest = classN; % the rest class is always the last

            trialsN = length(count_sequence(trial_channel,5));  % "5" tags each trial
            
        end;
        
        disp(['Classes/Movements detected: ' int2str(classN)]);
        disp(['Total trails detected: ' int2str(trialsN)]);
        disp(['Frequency: ' int2str(Fs)]);

        for m = 1:classN

            %classes should be supplied to this function
            %select only the trials for this class
            %"best_index" - in 1 shows the postion of a trial, "best_trigger" -
            %again shows the postion but in the range 1 .. n
            if (is_pulsar_single_trigger_channel_from_eeglab)
               %Pulsar GTEC
               [best_index,best_trigger] = generate_trials_postitions_using_single_discrete_channel(m, trigger_channel, 1, 2);
            else
               %Pulsar BIOPAC 
               [best_index,best_trigger] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            end;
            
            % cut signal
            %index = find(diff(trigger)==1); % because we do it earlier now!!!!
            trials_per_class = zeros(size(datar,2),window,length(best_index));
            disp(['Trials: ' int2str(length(best_index)) ' ' FileList{curFileIx}])

            %get the data for these triggers
            for i=1:length(best_index)
                %calculate a range around the trigger
                shift = 100;
                trials_per_class(:,:,i) = datar(best_index(i)+shift:best_index(i)+shift+window-1,:)';
            end

            X = cat(3,X,trials_per_class); % put all trials together
            Y = cat(1,Y,ones(length(best_index),1) * m);

        end;
   
    disp('--------------------------------');% next file
end

%% Classification over X and Y

% covariance estimation
COV = covariances(X);

%shuffuling
NTrial = size(X,3); % trials equals movements
ix = randperm(NTrial);

% permutate the data for classification
COV = COV(:,:,ix); % the sequence is updated from 1 ... ntrial to a random permutation
Ys = Y(ix); % it takes Y in the permutated form

% cross validation
NCV = 10;
% MDM
disp('---------- MDM --------------');
out = cross_valid(COV,Ys,NCV,@mdm); % COV is produced from X and Ys is produced from Y
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);


% FGMDM
disp('---------- FGMDM --------------');
out = cross_valid(COV,Ys,NCV,@fgmdm);
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);

% K-NN
disp('---------- KNN --------------');
% 10 neigbourg
out = cross_valid(COV,Ys,NCV,@knn,10);
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);

% LDA on std
disp('---------- LDA on std (no Riemannian geometry) --------------');
Feat = sqrt(squeeze(sum(X(:,:,ix).^2,2)))';
%  Constants
NTrial = size(X,3);
NTrialTest = fix(NTrial/NCV);
out = zeros(size(X,3),1);